#!/usr/bin/perl
#author: Bing Zhang

use strict;
use warnings;
use Statistics::Descriptive;
use Getopt::Std;

my %args;

my $usage = "pidMedian2eidStandard.pl -i <input gene expression matrix> -l <log2 transformation, 0/1, default 0> -m <pid 2 gid mapping file, \"NA\" if input id is gid> -t <mapping file type, wg, affy, soft, or ortholog> -p <pid column, count from 0, wg/0, soft/0, affy/0, ortholog/0> -e <eid column, count from 0, wg/2, ortholog/1, soft symbol/1, soft entrez 2,affy symbol/14, affy entrex/18> -s <standardize method: none, zscore, meandiff, mediandiff, or madz, default none> -c <columns used to calculate the means/medians and stdevs/mads for each gene, column numbers separated by \",\", optional, default all columns> -n <include columns used for standardization in the output, 0/1, default 0> -a <sample annotation file, optional> -o <output file>\n";

#file format:
#gene expression matrix file: output from rma etc, sample name column=sample data column-1
#pid2gid files: WebGestalt database table file, affymetrix annotation file, or ortholog mapped file (pid\teid\n)
#sample annotation file: sample by attribute matrix, sample name column=sample attribute column, sample name should match those in the gene expression matrix file

getopt("i:l:m:t:p:e:s:c:n:a:o", \%args);

my $expression_file = $args{i};
my $log2 = $args{l};
my $mapping_file = $args{m};
my $mapping_type = $args{t};
my $pid_column = $args{p};
my $eid_column = $args{e};
my $standardize = $args{s};
my $columns = $args{c};
my $include = $args{n};
my $ann_file=$args{a};
my $output_file= $args{o};

die $usage unless (defined ($expression_file) and defined ($mapping_file) and defined ($mapping_type) and defined ($pid_column) and defined ($eid_column) and defined ($output_file));

unless ($mapping_type eq 'wg' or $mapping_type eq 'ortholog' or $mapping_type eq 'soft' or $mapping_type eq 'affy'){
   print "Error: $mapping_type is not a valid mapping type. Supported mapping types are  wg , affy, and ortholog!\n";
   exit;
}

if (defined ($standardize)){
   if ($standardize ne "none" and $standardize ne "zscore" and $standardize ne "madz" and $standardize ne "mediandiff" and $standardize ne "meandiff"){
      print "Standardize method $standardize is not supported. Supported methods are none, zscore, meandiff, mediandiff, and madz.";
      exit;
   }
}
else{
   $standardize = "none";
}

$include = 0 unless (defined $include);
$log2=0 unless (defined $log2);

my %column_hash;
if (defined ($columns)){
   my @column_array=split(",",$columns);
   foreach my $column (@column_array){
      $column_hash{$column}=1;
   }
}

my %pid2eid;
my %ann2sid;

#read sample annotation files into a hash %ann2sid{ann}{sampleid}=annotationValue
if (defined ($ann_file)){
   my $linect=0;
   my @anns;
   open (ANNFILE, "$ann_file");
   while (my $line = <ANNFILE>){
      chomp $line;
      $linect++;
      if ($linect==1){
         @anns=split (/\t/, $line);
      }
      else{
         my @linearray = split (/\t/, $line);
         for (my $i=1; $i<@linearray; $i++){
            $ann2sid{$anns[$i-1]}{$linearray[0]}=$linearray[$i];
         }
      }
   }
   close (ANNFILE ); 
}     


if ($mapping_file ne "NA"){
#read pid2eid mapping file into a hash %pid2eid{pid}=eid
open (MAPPINGFILE, "$mapping_file");
while (my $line =<MAPPINGFILE>){
   chomp $line;
   my @linearray;
   if ($mapping_type eq 'wg'){
      @linearray=split (/\t/, $line);
      next if $linearray[1] >1;
      next unless $linearray[$eid_column];
      $pid2eid{$linearray[$pid_column]}=$linearray[$eid_column];
   }
   elsif ($mapping_type eq 'affy'){
      @linearray=split (/\"\,\"/, $line);
      next unless $linearray[$eid_column];
      $linearray[$pid_column] =~ s/^\"//;
      next if $linearray[$eid_column] eq '---';
      next if $linearray[$eid_column] =~ /\/\/\//;
      $pid2eid{$linearray[$pid_column]}=$linearray[$eid_column];
   }
   elsif ($mapping_type eq 'ortholog'){
      @linearray=split (/\t/, $line);
      next unless $linearray[$eid_column];
      $pid2eid{$linearray[$pid_column]}=$linearray[$eid_column];
   }
   elsif ($mapping_type eq 'soft'){
      @linearray=split (/\t/, $line);
      next unless $linearray[$eid_column];
      next if $linearray[$eid_column] eq '---';
      next if $linearray[$eid_column] =~ /\/\/\//;
      $pid2eid{$linearray[$pid_column]}=$linearray[$eid_column];
   }
}
close (MAPPINGFILE );
}

#read expression file into a hash or array @{%expression{eid}{sample}}=CorrespondingPidSignals
open (EFILE, "$expression_file");
my $i=1;
my @samples;
my %pidct;
my %expression;
while (my $line = <EFILE>){
   chomp $line;
   my @linearray = split (/\t/, $line);
   if ($i==1){
      $i++;
      for (my $j=0; $j<@linearray; $j++){
         my $sample = $linearray[$j];
         $sample =~ s/\.CEL//;
         push(@samples,$sample);
      }
   }
   else{
      if ($mapping_file ne "NA"){
         next unless exists $pid2eid{$linearray[0]};
      }
      for (my $j=1; $j<@linearray; $j++){
         $linearray[$j]=log($linearray[$j])/log(2) if ($log2==1);
         if ($mapping_file ne "NA"){
            push (@{$expression{$pid2eid{$linearray[0]}}{$samples[$j-1]}}, $linearray[$j]);
         }
         else{
            push (@{$expression{$linearray[0]}{$samples[$j-1]}}, $linearray[$j]);
         }
      }
   }
}
close (EFILE );


#unless (defined $columns){
#   for (my $l=1; $l<=@samples; $l++){
#      $column_hash{$l}=1;
#   }
#}

my $head;
open (OFILE, ">$output_file");
for (my $l=1; $l<=@samples; $l++){
   unless (exists $column_hash{$l} and $include==0){
      $head .= "$samples[$l-1]\t";
   }
}
$head =~ s/\t$//;
print (OFILE "$head\n");

if (defined ($ann_file)){
   foreach my $ann (sort keys %ann2sid){
      print (OFILE "sample_annotation_$ann");
      for (my $l=1; $l<=@samples; $l++){
         unless (exists $column_hash{$l} and $include==0){
            print (OFILE "\t$ann2sid{$ann}{$samples[$l-1]}");
         }
      }
      print (OFILE "\n");
   }
}

foreach my $eid (sort keys %expression){
   print (OFILE "$eid");
   my %expression_median;
   my $eid_signals = Statistics::Descriptive::Full->new();
   my @eidsignals;
   for (my $k=1; $k<=@samples; $k++){
      my $sample = $samples[$k-1]; 
      my $signals = Statistics::Descriptive::Full->new();
      foreach (@{$expression{$eid}{$sample}}){
         $signals->add_data($_);
      }
      my $median = $signals->median();
      $expression_median{$sample}=$median;
      if (!(defined $columns) or exists $column_hash{$k}){
         $eid_signals->add_data($median);
         push (@eidsignals, $median);
      }
   }
   if ($standardize eq "none"){
      for (my $l=1; $l<=@samples; $l++){
         unless (exists $column_hash{$l} and $include==0){
            my $ori_sig = $expression_median{$samples[$l-1]};
            $ori_sig = sprintf("%.4f",$ori_sig);
            print (OFILE "\t$ori_sig");
         }
      }
   }
   elsif ($standardize eq "zscore"){
      my $mean = $eid_signals->mean();
      my $stdev= $eid_signals->standard_deviation();
      $stdev=0.0000000001 if $stdev==0;
      for (my $l=1; $l<=@samples; $l++){
         unless (exists $column_hash{$l} and $include==0){
            my $zscore = ($expression_median{$samples[$l-1]}-$mean)/$stdev;
            $zscore = sprintf("%.4f",$zscore);
            print (OFILE "\t$zscore");
         }
      }
   }
   elsif ($standardize eq "meandiff"){
      my $mean = $eid_signals->mean();
      for (my $l=1; $l<=@samples; $l++){
         unless (exists $column_hash{$l} and $include==0){
            my $mean_diff = $expression_median{$samples[$l-1]}-$mean;
            $mean_diff = sprintf("%.4f",$mean_diff);
            print (OFILE "\t$mean_diff");
         }
      }
   }
   elsif ($standardize eq "mediandiff"){
      my $median = $eid_signals->median();
      for (my $l=1; $l<=@samples; $l++){
         unless (exists $column_hash{$l} and $include==0){
            my $median_diff = $expression_median{$samples[$l-1]}-$median;
            $median_diff = sprintf("%.4f",$median_diff);
            print (OFILE "\t$median_diff");
         }
      }
   }
   elsif ($standardize eq "madz"){
      my $median = $eid_signals->median();
      my $diff2median = Statistics::Descriptive::Full->new();
      foreach my $sig (@eidsignals){
         my $diff = abs($sig-$median);
         $diff2median->add_data($diff);
      }
      my $mad= $diff2median->median();
      for (my $l=1; $l<=@samples; $l++){
         unless (exists $column_hash{$l} and $include==0){
            my $madz = 0.6745*($expression_median{$samples[$l-1]}-$median)/$mad;
            $madz = sprintf("%.4f",$madz);
            print (OFILE "\t$madz");
         }
      }
   }
   print (OFILE "\n");
} 
close (OFILE );
